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ABSTRACT 

We present new dynamical models of the SO galaxy N3115, making use of the available 
published photometry and kinematics as well as of two-dimensional TIGER spectrog- 
raphy. 

The models are based on a detailed model of the luminosity distribution built using 
an MGE fit on HST/WFPC2 and ground-based photometric data. We first examined 
the kinematics in the central 40" in the light of two-integral f(E, J) models. Jeans 
equations were used to constrain the mass to light ratio, and the central dark mass 
whose existence was suggested by previous studies. The even part of the distribution 
function was then retrieved via the Hunter & Qian formalism. We thus confirmed 
that the velocity and dispersion profiles in the central region could be well fit with a 
two-integral model, given the presence of a central dark mass of ~ 10 9 Mq. However, 
no two-integral model could fit the /13 profile around a radius of about 25" where the 
outer disc dominates the surface brightness distribution. 

Three integral analytical models were therefore built using a Quadratic Program- 
ming technique. These models showed that three integral components do indeed pro- 
vide a reasonable fit to the kinematics, including the higher Gauss-Hermite moments. 
Again, models without a central dark mass failed to reproduce the observed kinemat- 
ics in the central arcseconds. This clearly supports the presence of a nuclear black hole 
of at least 6.5 x 10 8 Mqui the centre of NGC 3115. These models were finally used to 
estimate the importance of the dark matter in the outer part of NGC 3115, suggested 
by the flat stellar rotation curve observed by Capaccioli et al. (1993). 

This study finally points out the difficulty of integrating independently published 
data in a coherent and consistent way, thus demonstrating the importance of taking 
into account the details of the instrumental setup and the reduction processes. 
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1 INTRODUCTION 

Early-type disc galaxies pose a specific problem for dynam- 
ical modeling as neither the bulge nor the disc can be ne- 
glected in the computation of the gravitational potential. 
First, it is generally difficult to disentangle the photometric 
components according to a priori fixed criteria such as the el- 
lipticity or surface brightness profiles. Although there exists 
different methods to build corresponding realistic photomet- 
ric models, they are usually limited to two component sys- 
tems (e.g. disc plus spheroid, Byun & Freeman 1995) which 



* Based on observations taken with the Canada-France-Hawaii 
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do often not reflect the true morphology of the galaxy. Then, 
the different internal dynamics of each component has to 
be taken into account. Asymetrical Line Of Sight Velocity 
Distributions (LOSVDs) observed in discy ellipticals (e.g. 
Scorza & Bender 1995) confirmed the fact that the clas- 
sically measured velocity V and dispersion a do not suffi- 
ciently constrain their dynamical structure. Because of these 
difficulties, the analysis is usually restricted to the central 
region, or to regions where one component dominates the 
light (e.g. Jarvis & Freeman 1985). 

A number of techniques to build realistic distribution 
functions of such complex objects have recently been devel- 
oped. Hunter & Qian (1993) proposed a viable scheme to 
retrieve the even part of the axisymmetric two integral dis- 
tribution function corresponding to a given mass model. It 
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was used by van den Bosch & de Zeeuw (1996) for a set of 
models including a spheroid and a disc. In principle, this 
method can be applied to more complex potentials: this will 
be demonstrated in the present paper. As a further step, 
three integral axisymmetric models can be built, using for 
example the Schwarzschild method in which one derives the 
best combination of orbits constrained to fit a set of observa- 
tional data (e.g. van der Marel et al. 1998, Cretton & van den 
Bosch 1998). In this paper, we wish to present the applica- 
tion of another technique, namely Quadratic Programming 
(Dejonghe 1989) 

NGC 3115 is interesting in many different aspects. First 
it has long been assumed to be the prototype of the SO 
galaxy type: a bulge dominated galaxy with an embedded 
disc and very little gas and dust. With a distance of about 
10 Mpc and a nearly edge-on inclination, NGC 3115 is also 
a perfect target for optical observations. Numerous studies 
were focused on the understanding of its morphology and 
kinematics. NGC 3115 contains a double disc structure with 
an outer Freeman type II disc extending up to ~ 140" and a 
nuclear disc (inside 3") clearly revealed by the HST/WFPC2 
pictures (Kormendy et al. 1996). Its outer disc was found to 
exhibit a weak spiral arm structure (Capaccioli et al. 1987, 
Silva et al. 1988) and to thicken outwards (Capaccioli et 
al. 1988). The spheroidal component does not seem to follow 
the classical r 1 ^ 4 law, and the flattened halo of NGC 3115 
extends up to 20' (Capaccioli et al. 1987). 

The kinematics of NGC 3115 also deserved a complete 
survey: Illingworth & Schechter (1982) presented an exten- 
sive series of spectroscopic observations demonstrating that 
the bulge of this galaxy was nearly consistent with being 
an isotropic rotator. NGC 3115 has also been one of the 
best candidates for the presence of a central dark mass 
since the spectroscopic observations of Kormendy & Rich- 
stone (1992): spherical models suggested an increase of M/L 
by more than a factor of 10 inside 2", which is difficult to 
account for with normal stellar populations. This conclu- 
sion was based on the observed velocity gradient and cen- 
tral gaussian velocity dispersion ao, the latter nearly reach- 
ing 300 km.s _1 at a resolution of ~ 1" FWHM. This trend 
has been very recently confirmed with FOS spectroscopy 
(Kormendy et al. 1996), yielding a ~ 445 km.s _1 at HST 
resolution. Other studies mainly focused on the major- 
axis kinematics (Carter & Jenkins 1993, van der Marel et 
al. 1994, Bender et al. 1994, Fisher 1997). Finally Capac- 
cioli et al. (1993) have shown that the stellar rotation and 
dispersion profiles were flat at large radii (R > 100"), thus 
indicating the presence of a massive dark halo. 

In this paper, our aim is to make detailed models of this 
early type disc galaxy, e.g. constrain its distribution func- 
tion, fully taking into account the available published data. 
This includes, for the first time, two-dimensional spectro- 
scopic observations. In Sect. 2 we describe the photometry 
and spectroscopic observations used for the modeling, in- 
cluding new HRCAM images, the WFPC2 V band image 
as well as original TIGER spectrography of the central re- 
gion. The photometric models based on the MGE formalism 
are presented in Sect. 3 and their corresponding dynamical 
Jeans models in Sect. 4. Numerical two integral distribution 
functions are derived using the Hunter & Qian method and 
analysed in Sect. 5. Two and three integral models are corn- 



Table 1. Characteristics of the V band images of NGC 3115. 





Calar Alto 


HRCAM 


WFPC2 


pixel size 


1"55 


O'.'ll 


0'.'0455 


seeing (cr*) 


~3" 


0'.'7 




field of view 


10' x 15' 


70" x 110" 


34" x 33" 


exp. time 


240s 


120s 


1030s 


airmass 


1.65 


1.13 





puted in Sect. 6 with the help of the Quadratic Programming 
method. Conclusions are drawn in Sect. 7. 



2 OBSERVATIONAL DATA 

2.1 Photometry 

Dynamical modeling first implies the availability of a lu- 
minosity model. This requires images with high resolution 
to get details on the central region, but also with a large 
field of view to correctly sample the line-of-sight at large 
radii. A wide field V band image, obtained with a focal 
reducer on the 1.23m telescope at Calar Alto, was kindly 
put at our disposal by Cecilia Scorza. We then retrieved 
the V band WFPC2 image (F555W) of NGC 3115 from 
the HST archives at ECF (PI Faber). We finally used a 
V band image obtained by Jean-Luc Nieto in April 1992 
with HRCAM at the CFHT. All these images were reduced 
and normalized using aperture photometry available in the 
literature (Poulain 1986 & 1988, de Vaucouleurs & Longo 
1988). These calibrations led to individual errors smaller 
than 0.05 mag.arcsec~ 2 with an excellent overall agreement. 
The characteristics of the different images are given in Ta- 
ble §. 

2.2 Spectroscopy and kinematics 

2.2.1 TIGER data 

We observed NGC 3115 with the TIGER integral field spec- 
trograph (Bacon et al. 1995) at the CFHT in April 1992 and 
subsequently in April 1996. The lens diameter was set to 
0'.'39 and the spectral domain to 5130-5520 A for both runs. 
The spectra were extracted and calibrated in the standard 
way using the TIGER software written at the Lyon Obser- 
vatory. This includes bias subtraction, flat-fielding, spectra 
extraction, wavelength calibration, cosmic rays removal and 
differential refraction correction (see Emsellem et al. 1996 
for more details). The individual exposures were then accu- 
rately centred and merged. Variations in the overall atmo- 
spheric transparency were smaller than 1% and corrected. 
We thus obtained one set of spectra for each run, consisting 
of 624 and 391 spectra for the run92 and run96 respectively. 
The final seeing of each data cube was estimated a poste- 
riori by comparing images reconstructed from the spectra 
with the V band WFPC2 image. Details on the character- 
istics of both sets of data are given in Table |^. Note that 
although the set of data obtained in April 1996 has a signifi- 
cantly better spatial and spectral resolution, we will mainly 
present the data set obtained in April 1992 as the signal to 
noise is higher and the field of view larger. 
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Figure 1. TIGER stellar kinematical maps: velocity (top left), dispersion (top right), /13 (bottom left) and /14 (bottom right). The 
isocontour step is 25 km.s" 1 for both V and <r, and 0.05 for /13 and hi. 



Table 2. Observational characteristics of the TIGER spectro- 
graphic exposures. The fields of view are only indicative as the 
coverage of the merged TIGER spectra are not rectangular. All 
individual frames had an exposure time of 30 mn. 





1992 Run 


1996 Run 


Lens diam. 


0^39 


0'.'39 


Field of view 


14'.'5 x 7'.'9 


7'.' 9 x 7('9 


# of exposures 


4 


2 


# of spectra 


624 


391 


Spect. sampling 


1.8 A .pixel -1 


1.5 A .pixel" 1 


Spect. resolution 


3.7 A (FWHM) 


3.2 A (FWHM) 


Instr. broadening (<r) 


94 km.s" 1 


78 km.s" 1 


Spect. domain 


5130 - 5520 A 


5130 - 5520 A 


Seeing (FWHM) 


T/24 


0'.'8 


Stellar templates 


7] Cygni 


HR 127065 




K0III 


K0III 



In order to determine the stellar kinematics we also ob- 
tained exposures of a few stellar templates: mainly K giants 
for run92 and G, K and M giants for run96. All spectra 
were logarithmically rebinned, the stellar continuum was ap- 
proximated using a fifth degree polynomial and subtracted. 
We then derived the 2D kinematics of the central region 



of NGC 3115 using the FCQ algorithm (Bender 1990) and 
an optimal template built with the available stellar spectra. 
This provides the full LOSVD for each spatial element. The 
LOSVDs were then parametrized with the Gauss-Hermite 
functions. True moments were also estimated from the pos- 
itive part of the parametrized LOSVDs (see Emsellem et 
al. 1996). The results were found to be rather insensitive to 
the choice of the template, so we decided to use the spec- 
trum of a single K0 giant for both runs. We finally built the 
maps for each kinematical quantity: V , a, /13 and hi which 
correspond to the Gauss-Hermite parametrization as well as 
the estimations of the four first true moments of the LOSVD 
V , a, the skewness and the kurtosis. 

The mean velocity, velocity dispersion (gaussian fit 
/13 and hi two-dimensional maps are presented in Fig. 
(run92). We detect a slight tilt of 10° of the velocity minor- 
axis in the 1992 data. The low signal to noise of the new 
TIGER data set prevented us to confirm this: it would there- 
fore be useful to obtain new high spatial resolution kinemat- 
ical data along the minor-axis. Note that the /13 and hi maps 
only include spectra with a signal to noise higher than 20. 
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Figure 2. Comparison between the different published kinematics and the TIGER 1992 data. Each symbol corresponds to a reference 
as shown in the V panel. F97 = Fisher 1997; vdM = van der Marel et al. 1994 (WHT data only); BSG = Bender et al. 1994; KR = 
Kormendy & Richstone 1992 (exposures 29F48, 29F42, 34F16; see KR92); CJ = Carter & Jenkins 1993. 



Table 3. Observational characteristics of the kinematical (long-slit) data used in this paper (for the FOS/HST data, see Kormendy et 
al. 1996). 



Ref. 



slit width 



seeing (FWHM) h 3 ,h 4 



Illingworth & Schechter 1982 


2-3" 


100-130 kirns" 1 


l'.'5-5" 


no 


Capaccioli et al. 1993 


l-T/6 


~ 40 km.s" 1 


> l'/5 


no 


Kormendy & Richstone 1992 


(f.' 5 


68 km.s -1 


1"-T.'6 


no 


Carter & Jenkins 1993 


C'45 


10 krn.s" 1 


0'.'8 


no 


van der Marel et al. 1994 


1" 


38 km.s -1 


0'.'96 


yes 


Bender et al. 1994 


2'.'1 


46 km.s -1 


1-2'.' 5 


yes 


Fisher 1997 


2" 


75 km.s -1 


> 1'.'5 


yes 


Kormendy et al. 1996 


r/.'26 


~ 35 km.s 


0'.'57 





2.2.2 Comparison with published kinematics 

All kinematical data available were compiled, and the 
FOS/HST (0'.'21 square aperture) and SIS/CFHT data ob- 
tained by Kormendy et al. (1996) scanned. As emphasized 
by different authors (e.g. van der Marel et al. 1994) the gaus- 



sian velocity V and dispersion a are not always good esti- 
mations of the true V and a moments. This is particularly 
important in the case of NGC 3115 where high values of /13 
have been measured (e.g. van der Marel et al. 1994). There- 
fore, we mainly focused on the data sets including higher 
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order moments of the LOSVDs: van der Marel et al. (1994, 
hereafter vdM+94), Bender et al. (1994, hereafter B+94), 
Fisher (1997, hereafter F97) and the TIGER data. In these 
cases, we used the Gauss-Hermite expansion of the LOSVDs 
to estimate V and a using the same method as used by Em- 
sellem et al. (1996). Still we made use of other published 
kinematics such as Kormendy & Richstone (1992, hereafter 
KR92), Carter & Jenkins (1993), and the recent high spa- 
tial resolution FOS/HST and SIS/CFHT data (Kormendy 
et al. 1996, hereafter K+96). We finally included the data 
of Illingworth & Schechter (1982, hereafter IS82) for offset 
axes and of Capaccioli et al. (1993, hereafter C+93) for the 
kinematics at large radii along the major-axis. The charac- 
teristics of each data set are given in Table ^. 

In Fig. ^, we compare the TIGER kinematics along the 
major-axis with the ground-based data with similar resolu- 
tion: they are all in good agreement. The error bars of the 
TIGER data were determined by the statistics on the fluctu- 
ations inside a spatial resolution element, here taken as the 
FWHM of the PSFp|. As it turns out, this defines a beam 
containing on average 7 spectra. 



2.3 Convolution and pixel integration 

The comparison between observed kinematics and theoret- 
ical models requires to take into account the observational 
characteristics. The seeing smearing and pixel binning have 
been computed through a double quadrature (see Appendix 
A). All the comparisons presented in the following sections 
include this processing, according to the parameters of the 
relevant data set as tabulated in Table ^. Note that the 
PSF of the FOS data was derived using an approximation 
as given in van der Marel et al. (1997). 

2.4 Uncertainties in the kinematics 

Additional factors linked with the instrumental setup or 
the data reduction processes can significantly influence the 
measured kinematics, and in particular the higher order 
Gauss-Hermite moments. Firstly, an error in the center- 
ing or position angle of the slit induces a change in the 
observed kinematics. In the case of NGC 3115, an offset 
of Aa = 2° (major-axis) in the position angle of the slit 
gives A/13 ~ 10% at R = 35", while the change in ve- 
locity amounts to only about 3%. It is also important to 
control the spectral filtering applied for the retrieval of the 
LOSVDs as it could lead to underestimate the higher order 
moments. Finally, the obtained kinematics may be sensitive 
to the method used (e.g. Fourier fitting, FCQ, etc). In the 
following paragraph, we discuss such an effect on the central 
FOS LOSVD. 



2.4-1 The central hi value 

The central observed FOS LOSVD is one of the best data set 
we can use to constrain the central dark mass in NGC 3115, 

T Lens to lens fluctuations in the reconstructed kinematical maps 
are a good estimator of the errors, formal and instrumental 



as it has the highest available spatial resolution, it in- 
cludes the higher order kinematical terms and does not de- 
pend on the odd part of the distribution function (or very 
weakly since r — O'.'Ol). K+96 quotes a central value of 
hi — 0.10 ±0.03. However, higher order Gauss-Hermite mo- 
ments may significantly depend on the instrumental setup 
and reduction procedures. We thus wish here to assess the 
robustness of the central hi value before attempting any 
detailed modeling. 

We have therefore retrieved and reduced the FOS 
spectra and tried to reproduce the analysis mentioned in 
K+96. Our velocity and dispersion profiles obtained with 
the Fourier Quotient method (FQ) are compatible (within 
15 km.s -1 ) with the ones quoted by K+96: e.g. we find a 
value of cr(0) = 448 km.s -1 to be compared with cr(0) = 
443+18 km.s -1 as given by K+96. However, the shape of the 
LOSVD, and thus the dispersion and hi values, obtained us- 
ing the Fourier Correlation Quotient (FCQ) can be severely 
influenced by template mismatch, a bad continuum subtrac- 
tion and the spectral filtering (see also vdM+94). These 
would in turn affect the high frequencies of the LOSVD, 
thus diminishing its peakedness, and/or the low frequencies 
smoothing out its wings. Moreover, large hi values may be 
anyway difficult to obtain, as already emphasized by van der 
Marel (1994). 

We have thus simulated artificial FOS spectra by con- 
volving a FOS stellar template with different numerical 
LOSVDs, including the noise level present in the observed 
FOS spectraQof NGC 3115. We then calculated the LOSVDs 
via FCQ, as in K+96, and parametrized them using Gauss- 
Hermite functions. These simulations clearly demonstrate 
that large hi values are underestimated (and as a con- 
sequence, gaussian dispersion obtained via FCQ proba- 
bly overestimated): the dominant effect is a strong filter- 
ing of the central peak, although the large wings of the 
LOSVDs are also significantly affected. Considering the cen- 
tral LOSVD presented in K+96 and the results of our sim- 
ulations, we can estimate the true value of the central hi 
to be ~ 0.17 with a possible range 0.14 < h 4 < 0.21. But 
as discussed above, these values are ill-constrained and an 
accurate estimation would require additional observations. 
We should finally note that this sensitivity is particularly 
important for the FOS/HST data: at a lower spatial resolu- 
tion, gradients are not as extreme and the LOSVDs closer 
to Gaussians. 



3 THE MGE LUMINOSITY MODEL 

We applied the MGE technique (Monnet et al. 1992, Em- 
sellem et al. 1994, Emsellem 1995) to the available V band 
images to build a complete photometric model of NGC 3115. 
The individual Point Spread Functions were taken into ac- 
count to obtain a deconvolved model which fits the data 
from the central HST/WFPC2 pixel (OC'0455) to 300". In 
Fig. ^ we present the isophotes of NGC 3115 at three differ- 
ent scales compared with the MGE V band model. The fit 
is excellent at all scales (see also Fig.^). The only significant 

t This noise does not appear in Fig. 3 of K+96, so we assumed 
that the spectra have been significantly filtered for the plot. 
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Figure 3. Isophotes of the V band images of NGC 3115 
(thin lines) and of the corresponding convolved MGE model 
(19 components): wide field image (Calar Alto, top), HRCAM 
(middle) and WFPC2 (bottom). The isophotes step is 
0.5 mag.arcsec - 2 , and the faintest isophotes are 25, 20.5 and 
17.5 mag.arcsec - 2 respectively 



discrepancies are observed at R ~ 5" where the isophotes are 
boxy due to a detected spiral-like structure, and at R ~ 30" 
where the galaxy slightly departs from axisymmetry. The 
largest component has a FWHM of 895" and an axis ratio 
of 0.85. 

We did not attempt to go beyond this scale as the photom- 
etry becomes rather uncertain. Hence, at radii larger than 
600" along the major-axis, the MGE model decreases more 
rapidly (exponentially) than the B band profile published by 
Capaccioli et al. (1987). However we tested that this does 
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Figure 4. Central V band major and minor-axis profiles. The 
HST/WFPC2 points are compared with the MGE models and a 
three integral model (dashed curves, see Sect. along the major- 
axis (left panel) and minor-axis (right panel). In each panel, the 
lower solid and dotted curves correspond to the convolved (and 
binned) MGE models with a central gaussian (solid lines) or a 
cusp (7 = 1.5, dotted lines; see Appendix B). The corresponding 
deconvolved MGE models are also shown for comparison (two 
upper curves in each panel). 



not significantly influence the dynamics up to 150", radius 
of the last observed kinematic data point. 

In Fig. we show the isophotes of the deprojected MGE 
model for an inclination angle of 86°. The different compo- 
nents clearly appear in this plot: 

• the flattened extended halo which dominates the light 
of the galaxy at radii larger than 100" (e ~ 0.65 at 300"). 

• the outer disc which extends up to ~ 100" and is trun- 
cated at R ~ 5" (Freeman type II disc); 

• the inner spheroid which exhibits a boxy shape and an 
ellipticity in the range 0.4 — 0.6; 

• the inner or nuclear disc which continues up the the 
centre and has a radial extent of about 3"; 

• and finally the point-like nuclear source with a mag- 
nitude mv = 16.63 mag.arcsec -2 and a half light radius of 
vh ~ 0"054, values consistent with the ones given by Kor- 
mendy et al. (1996). Note that it is yet impossible to disen- 
tangle between a flattening of the surface brigtness profile 
of this nucleus inside 0.1" and a still increasing power law 
(Fig.§. 

In the following modeling we will always use a distance of 
10 Mpc for NGC 3115: this value is intermediate between 
11.0 recently given by Elson (1997) and the previously used 
one (9.24 Mpc - e.g. KR92). 



4 TWO INTEGRAL JEANS MODELS 

In this Section we wish to constrain the mass to light ratio 
in the central 40" of NGC 3115 by comparing the observed 
kinematics with Jeans dynamical models. We assumed a 
constant M/Ly for the luminous mass of the galaxy, and 
added any required additional dark mass. We computed 
the second order non-centred projected velocity moment 

[12 = \/v 2 + <t 2 (V and a are the projected mean velocity 
and velocity dispersion respectively) by solving the Jeans 
equations for a two-integral distribution function f(E,J). 
This is particularly straightforward in the context of the 
MGE formalism as /12 can be evaluated through a single 
quadrature (see Emsellem et al. 1994). Moreover, in the 
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Figure 5. Isophotes of the deprojected (and deconvolved) MGE 
model of NGC 3115 (V band) in the meridian plane (R, z): the 
steps are 0.5 in log (ijQ.pc -3 ), with the faintest isophotes corre- 
sponding to (from top to bottom) -4, -1.5, and 0.5 respectively. 



The nucleus peaks at 
central gaussian). 



5.42 (~ 1.29 10' I^.pc -3 , model with a 



frame of a f(E, J) dynamical model, \i2 does not depend 
on the dispersion tensor anisotropy. Observed [12 can only 
be determined with reasonable accuracy for data including 
higher order moments. 

4.1 The inclination angle 

Taking into account the axis ratios of the disc components, 
the inclination angle i of the galaxy is restricted to a small 
range between about 83° and 90° (edge-on). Note that the 
lower limit is slightly model dependent as the MGE model 
cannot exactly reproduce the thickness of the outer disc 
which is constant or slightly rising outwards (see Capaccioli 
et al. 1988). We have built Jeans models assuming differ- 
ent values for i, sampling the interval between 83° and 90°. 
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Figure 6. Major and minor axes second order moment profile 
(squares: vdM+94; circles: F+96; losanges: KR92; crosses: IS82) 
compared with the MGE jeans model with M/Ly = 6.5 (solid 
lines) and M/Ly = 5.8. 



Varying the value of i mostly affects the major- axis kine- 
matics, where the flattened components dominate the light 
distribution. The minor-axis \i2 profile remains nearly un- 
changed. We can thus use the observed ratio between the 
value of /i2 along the minor-axis and along the major-axis 
as an indicator for the inclination angle, as it does not de- 
pend on the mass to light ratio and anisotropy. A value of 
i = 86° gives the best fit, and will be used for all following 
models. 



4.2 The stellar mass to light ratio 

Fig. |^ clearly shows that the central second order moment 
is underestimated in the model with constant M/Ly. In 
the restrictive frame of our axisymmetric Jeans models with 
constant M/L that we consider in this section, we could not 
find a reasonable fit to the data in the central 3". Outwards 
from the central 3", we obtained a rather good agreement 
with the [12 profiles along both the minor and major axes 
using M/Ly = 6.5 (Fig. 0). 



4.3 An estimation of the central dark mass 

As emphasized already in the previous Section and by pre- 
vious studies (see K+96 and references therein), the central 
dispersion gradient calls for an increase of the mass to light 
ratio in the central few arcseconds of NGC 3115. The FOS 
kinematics presented by K+96 then implies that this takes 
place inside a radius of ~ 0'.'2. According to our Jeans mod- 
els, the observed gradients in the kinematics would lead to 
M/Ly greater than 60 for the nucleus. This strongly sug- 
gests the presence of a central dark mass, which we denote 
by Mbh, that is larger than 10 8 Mq (the mass of the nu- 
cleus with M/Ly = 6.5 would be ~ 6.3 x 10 7 Mq). In the 
following paragraphs we now provide some constraints on 
the central mass concentration in the frame of a f(E, J) 
dynamical model. 

The addition of a central point-like mass in the MGE 
model is straightforward and the formalism is described in 
Appendix A of Emsellem et al. (1994). We built a number 
of models with masses Mbh ranging from 10 s to 2 x 10 9 Mq 
and for different M/Ly . In the following, we give the values 
of A^bh which best fit each individual data set, as well as 
an indicative range in an attempt to take into account the 
formal error bars of the measured kinematics. Note that fit- 
ting the central kinematics requires to vary both Mbh and 
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Figure 7. Observed kinematics for 3 data sets (upper panel: 
vdM+94, middle panel: TIGER, lower panel: B+94) along the 
major-axis compared with models with different black hole masses 
A/bh: 0, 3.7 x 10 8 , 6.5 x 10 s , 1.0 x 10 9 and 2.0 x 10 9 Mq. 



M/Lv ■ Fig. presents the comparison between these models 
and data inside 7" along the major-axis. 

Estimations of the true moments are sensitive to errors 
in the measured high order Gauss-Hermite moments (e.g. 
/13 and hi, see Bender et al. 1994): this may be the cause 
of the low central dispersion value in vdM+94's data. If we 
ignore the two central points in the fit we find Mbh ~ 6 ± 
1 x 10 8 Mq and M/Lv = 6.1 ± 0.4. The overall best fit to 
the vdM+94 data is obtained with M BH ~ 4± 1.5 x 10 8 M 
and M/L v = 6.2T0.3. The TIGER data indicates a slightly 
higher value for the central dark mass of A/bh = 6.6 ± 1 x 
10 8 Mq with M/Lv = 6.5 ± 0.5. Similar comparisons have 
been done with the data of Fisher (1997) and Bender et 
al. (1994) which both give similar values of Mbh ~ 7.0 ±3 x 
10 8 Mq and M/L v = 6.8 ± 0.5. 



For the global set of kinematical data, the best estima- 
tion for the black hole mass is then A/bh = 6.5±3.5x 10 s Mq 
with M/Lv = 6.5 ±0.7. Values as high as 2 x 10 9 quoted by 
K+96 are excluded at more than the 3<r level for the simple 
f(E, J) casd 



4.4 Dark matter at large radii 

At radii larger than 40" along the major-axis the observed 
velocity and dispersion profiles are almost flat out to the 
last measured point with (V) ~ 250 km.s _1 and (<r) ~ 
100 km.s -1 (Capaccioli et al. 1993). The Jeans model with 
M/Lv = 6.5 predicts a dispersion profile that is consistent 
with the observed one at large radii, but underestimates the 
velocity profile for R > 70" by as much as 60% at R = 120". 
As will be seen in Section [| this implies a strong increase 
of the mass to light ratio at large radii. 



TWO INTEGRAL DISTRIBUTION 
FUNCTIONS 



5.1 The method 



In Sect. ||, we have restricted the range of values for the 
central dark mass using Jeans models and estimations of 
the true moments, and found a best fit with Mbh = 6.5 ± 
3.5 x 10 8 Mq. We wish now to compare the models directly 
to the measured values, namely V , a and the higher order 
Gauss-Hermite moments. This requires the knowledge of the 
full distribution function. We therefore applied the method 
of Hunter & Qian (1993) to derive the even part of the distri- 
bution function f(E, J) corresponding to the axisymmetric 
MGE mass model of NGC 3115. 

The even part f e of f(E, J) was calculated with a grid 
of 49 x 38 points in (E, J)-space designed to properly sam- 
ple the gradients of / (particularly at J '(E) / 'Jmax(E) ~ 1 
due to the presence of the discs). The odd part f remains 
as a free function in our models with the constraint that 
f(E,J) = f e (E,J) + f (E, J) is positive everywhere. We 
have separated the contributions of the discs and the bulge in 
the computation of the distribution function: this is achieved 
by using the mass density of each component but including 
the total potential. 

In order to derive the odd part f a we used the 
parametrization introduced by van der Marel et al. (1994) 
following the work of Dejonghe (1986,1987): we define a 
function h a (r] = J / Jmax(E)) such as 



r tanh(a7?/2)/tanh(a/2) (a > 0) 

h a (ri) = \ n (o = 0) 

t (2/a) arctanh(»jtanh(a/2)) (a < 0) 



(1) 



so that f (E,J) — h a (rf)f e (E, J). In practice, we used dif- 
ferent values of a for the disc and bulge components. This 
was essential to fit the detailed observed kinematics. 

The LOSVDs were derived from the obtained numerical 
distribution function on a fine grid of more than 1000 points 
to allow a proper sampling of the sky plane in the central 



3 The distance used in this paper for NGC 3115 is slightly higher 
than the one of K+96, which increases this discrepancy. 
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Figure 8. Even part of the distribution function obtained with 
a mass model of NGC 3115 including a central dark mass of 
Mbh = 1-0 10 9 Mq: the logarithm of the distribution function 
is plotted as a function of the normalized energy E/E m ax and 
angular momentum J/Jmax- E max corresponds to the central 
potential of the model when the central dark mass is excluded. 
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Figure 9. Observed kinematics (losanges) along the major-axis 
(two top panels) and the minor-axis (bottom panels) for the data 
of KR92 and the best fit model with M BH = 0.94 X 10 9 Mq. 



45". The areas of the projected LOSVDs (normalized by the 
projected luminosity density) were always equal to 1 within 
1%: this is an a posteriori check of the validity of the over- 
all computations. This grid of LOSVDs was then used to 
compute (via an interpolation) the corresponding LOSVDs 
for each data set including the effects of seeing convolution, 
pixel integration and spectral resolution. Finally the result- 
ing LOSVDs were parametrized using the Gauss-Hermite 
moments as well as the true velocity moments. 

A number of models have thus been built, including dif- 
ferent cusp slopes and a range of values for the mass to light 
ratio M/Lv and for the central dark mass Mbh. In Fig. [s| 
we present one example of such a distribution function ob- 
tained for p oc r -1,5 and A/bh = 1.0 x 10 9 MQas a function 
of the normalized energy (E max corresponds to the central 
potential of the model excluding the central dark mass) and 
angular momentum (J = J ma x corresponds to circular or- 
bits at a given energy). The double disc structure is clearly 
visible as two bumps near the line of maximum angular mo- 
mentum, and the plateau at high energy corresponds to the 
cusp (for a cusp slope 7 = —1.5 and a keplerian potential 
the distribution function is constant with energy, see Qian 
et al. 1995). 



5.2 The inner 45 arcseconds 

As shown in Sect. 5|, an overall good fit to the data inside 40" 
requires a rather high mass to light ratio M/Lv ~ 6.5. For 
the best sets of data (namely vdM+94, CJ93, and KR92) 
which cover these radii, we built models with different cen- 
tral dark mass concentrations. In Fig. ^| we present the 
KR92 data set which extends further than 40" and over- 
imposed the best fit model which includes a central dark 



mass of Mbh = 0.94 x 10 9 MQand has a mass to light ratio 
M/Lv ~ 6.1. The major-axis velocity and dispersion (best 
gaussian V and a) profiles are well fit by this model even 
in the inner part. The fit is also good along the minor-axis 
besides a significantly lower dispersion around 3" where the 
central dark mass starts to dominate the observed kinemat- 
ics. The fit to the data of CJ93, which has a slightly higher 
spatial resolution, is also excellent. 

In Figure [ji^ we present the fit to the vdM+94 data set. 
In this case, the higher order Gauss-Hermite moments could 
also be included. We could not obtain a reasonable fit to the 
hz parameter for which the 21 model predicts values around 
-0.25 at R ~ 25" while the observed one is ~ —0.16. It is 
possible to reduce (in absolute value) the predicted /13 by 
changing the odd part of the distribution function, but at 
the expense of a significantly lower velocity V. Note that 
the predicted /13 does not depend on the mass to light ratio 
(assuming it is constant with radius). 

Uncertainties due to errors in the positioning of the slit 
or to spectral filtering (see Sect. 2.4) of the LOSVDs are 



far too small to fully account for this discrepancy. The ef- 
fect of the latter is shown in Fig. |l(j where the predicted 
LOSVDs have been slightly Fourier filtered before measur- 
ing the Gauss-Hermite moments: the fit is then significantly 
better for all /13, /14, /15, he profiles, without changing V and 
a significantly. However, the filtering required to explain the 
data is large considering the signal to noise (S/N ~ 20 per 10 
km.s -1 ) and spectral resolution (<Tmst ~ 23 km.s -1 ) quoted 
by vdM+94. We then performed some extensive simulations 
to derive the effect of the deconvolution method (Fourier 
fitting of the LOSVD) on V and I13: it is indeed significant 
for R ~ 25", where it can reach 15 km.s _1 and ~ 0.04 re- 
spectively. The data of vdM+94 are thus consistent with 
/13 = —0.2 at R ~ 25", still too low in absolute value to be 
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Figure 10. Observed kinematics (open squares) along the major- 
axis for the data published in vdM+94, derived using the method 
described by Rix & White (1993). These data have been fit simul- 
taneously on both sides of the centre (see vdM+94 for details). 
The best fit 21 model with M B h = 0.94 X 10 9 Mq is superimposed 
with (dotted line) and without (solid line) spectral filtering. 
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Figure 11. Observed kinematics (open squares) along the major- 
axis for the data of vdM+94, derived from a Fourier fitting 
method. Here the data has been derived on each side of the centre 
independently (north-east side: open triangles; south-west side: 
open squares). The best fit 21 model with M B h = 0.94 X 10 9 Mq 
is superimposed (solid lines). 



reconciled with our two-integral model. There are finally un- 
certainties associated with the data as illustrated in Fig. 
the same 21 model is compared with the kinematics obtained 
by vdM+94 but this time derived using a one side Fourier- 
fitting technique (see vdM+94 for details). The fits to the 
higher order moments are good, although the predicted V/cr 
is now too low. 



5.3 The central arcseconds 

The observed kinematics in the central arcseconds is very 
sensitive to the spatial resolution, particularly in the case 
of NGC 3115 for which both the photometry and kinemat- 
ics ex hibit large central gradients. But as we have shown in 
Sect. 2.4.1, the procedures used in the analysis of the data 



do also affect the final result. In the following paragraphs, 
we used the data provided by the FOS/HST and SIS spec- 
trographs in parallel with the (lower resolution but) two- 
dimensional TIGER kinematics to constrain the free param- 
eters of the dynamical models. We again examined a large 
range of different models, varying the central dark mass, the 
mass to light ratio, the central cusp slope as well as the odd 
part of the two- integral distribution function. For reasons of 
clarity, we only present here a few of these models which fit 
the data well. 



5. 3. 1 Best fit models 

All our two integral models predict rather high hi values 
(> 0.08) when M BH is a few 10 s Mq. Our best fit model 
to the FOS data only has a cusp with p oc r -1 ' 5 and a flat- 
tening of 0.7, a central mass M BH = 1.0 x 10 9 Mq and 
M/L v = 6.5 (model Ml, see Fig. [u|, although the cusp 



slope and flattening are not well constrained by the present 
data. In Fig. [l3], we present the kinematical profiles obtained 
from the FOS spectra and model Ml. The overall agree- 
ment is good, but our two integral model predicts a slightly 
higher central velocity dispersion with <r(0) = 480 km.s. . 
We should note however that the best gaussian fit to the 
LOSVD obtained by K+96, and presented in their Fig. 3 
has cr(0) ~ 489 km.s -1 , significantly larger than the one de- 
rived via FQ, but consistent with our value. This point is 
emphasized in Fig. ^ where we now present the comparison 
between the FOS LOSVD obtained by K+96 and the one 
predicted by model Ml. The agreement is excellent, consid- 
ering that the wiggles present in the FOS central LOSVD 
are almost certainly not real (as are the negative points), 
but resulting from the presence of noise in the FOS spec- 
trum as well as the template mismatching. It is striking to 
see the difference between the LOSVD directly predicted 
from model Ml and the one retrieved via FCQ using the 
same spectral filtering than for the FOS data: this forces us 
here to underestimate hi and correspo nding ly to overesti- 
mate the gaussian dispersion (see Sect. 2.4.1). 



Model Ml inc lud es a rapidly rotating nuclear disc with 
a = 5 (see Section 5.1). This model predicts velocities which 
are clearly too high at the SIS resolution, as shown in Fig. |l4| 
Kormendy et al. (1998) noticed the fact that, for SIS data of 
NGC 3377, FCQ seems to provide slightly higher velocities 
than FQ by about 10% (e.g. see their Fig. 15). If this is also 
the case for NGC 3115, it could reconcile our model Ml with 
the SIS data. 

We calculated other models where the parameter a was 
significantly increased (up to a = 100) for stars with high 
energy (close to the centre as we deal with positive energies) : 
in this region the total distibution function tends towards 
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Figure 12. Central LOSVDs: observed FOS (thin solid line, 
shifted to zero velocity), directly predicted from model Ml (dot- 
ted line), and the same obtained via FCQ and using the charac- 
teristics of the FOS central spectrum (thick solid line). 
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Figure 13. Observed FOS kinematics along the major-axis 
(K+96, crosses), compared with the velocity and dispersion pro- 
files predicted by models Ml, M2 and M3 (thin solid, dotted and 
thick solid lines respectively) via FCQ and using the characteris- 
tics of the FOS data. 
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Figure 14. Observed SIS kinematics along the major-axis 
(K+96, circles), with the thin solid, dashed and thick solid lines 
corresponding to predictions of models Ml, M2 and M3 respec- 
tively. 
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Figure 15. Predicted 2D kinematics for the TIGER, setup from 
model Ml. From top to bottom, left to right: V, cr, /13 and /14. 
Contours are the same than in Fig. l| 



the maximum streaming two-integral model. Fig. [13| and 
Fig. |lj present two other models, respectively M2 which cor- 
responds to Mbh = 1.24 10 9 M Q and M/L v = 5.4, and M3 
which has Mbh = 0.94 10 9 M and M/L v = 6.1. Although 
model M2 is a better fit to the FOS velocity and disper- 
sion profiles, the higher black hole mass implies hi = 0.134, 
a value at the upper limit of the error bar mentioned by 
K+96. Model M3 is the best fit to the SIS data, and repre- 
sents the best compromise to the HST data with a measured 
central /14 = 0.094. These models also show that the detailed 
kinematical structure inside 0'.'3 remains ill constrained. Al- 
though stellar orbits may be biased towards near maximum 
angular momentum values in the central 10 pc, it is thus yet 
improper to discuss the internal kinematical structure of the 
point-like nucleus detected in the photometry (Sec. ^|), all 
the more since its luminosity distribution is uncertain as well 
(see Fig. H). 



5.3.2 The TIGER data 

We finally compare the three best fits (Ml, M2 and M3) to 
the kinematics obtained with TIGER. Although the resolu- 
tion of these data is significantly lower than the SIS data, it 
still represents a good constraint as it homogeneously cov- 
ers a two-dimensional field on the sky. As seen in Fig. [j], the 
TIGER isovelocities are significantly flattened. This suggests 
again that the nuclear disc is rapidly rotating. The con- 
trast between the disc and the central spheroid is however 
weaker than in the case of the Sombrero galaxy for which 
very asymmetric LOSVDs (and therefore high /13 values) are 
observed at the centre (Wagner et al. 1989; Emsellem et al. 
1996). All three models, Ml, M2 and M3, fit reasonably well 
the TIGER kinematical maps, with Ml giving the smallest 
residuals. The kinematical maps corresponding to model Ml 
are presented in Fig. [l^. 

The results presented in the previous Sections, using 
two integral f(E, J) models to fit the observed kinematics, 
can thus be summarized in a few points: 

• There is a significant discrepancy between the observed 
and predicted /13 values at a radius of about 25". 
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• The fit of the central FOS kinematics requires the nu- 
clear disc to have a nearly maximum streaming distribution 
function inside l". 

• All the available kinematics in the central arcseconds 
can be rather well fit by a two-integral model with constant 
mass to light ratio. The best constrained parameters are the 
central dark mass and the mass to light ratio: considering 
the error bars in the different data sets, our best estima- 
tions are Mbh = 0.94+^j 5 10 9 M and M/L v = 6.1±£f. 
These values can be refined using the higher order Gauss 
Hermite moments in the central arcsecond. However the ob- 
served kinematics are very sensitive to the details of the 
instrumental setup and the analysis procedures. 

We could in principle reconcile our 21 model with the 
observed kinematics along the major-axis. The input mass 
model has been indeed derived from a V band image which 
may not reflect the true mass distribution. We have for ex- 
ample assumed that the bulge and disc components have the 
same mass to light ratio. This is certainly not true since they 
do not share the same colour profiles (Fisher et al. 1996). A 
difference of 0.1 magnitude in B — R c as typically observed 
between the bulge and the disc regions in NGC 3115 (Fisher 
et al. 1996) corresponds to a maximum change of about 5% 
in M/Lv (Worthey et al. 1994), too small an effect to ac- 
count for the discrepancy in hg. In the following, we wish 
to examine another explanation: the distribution function 
depends on a third integral of motion. A two integral distri- 
bution function f(E, J) imposes that or = a z everywhere. 
This may not be the case for the SO galaxy NGC 3115, as 
for our own Galaxy. This is also supported by derivations 
kindly made by Roeland van der Marel, using the formalism 
described in de Bruijne et al. (1996) which allows a rather 
general geometry of the velocity ellipsoid. Although these 
models cannot mimic the complex morphology of NGC 3115, 
they clearly suggest that, for very flattened components, low 
/13 values of the order of —0.2 could be obtained by relaxing 
the constraint of a two-integral distribution function. 



6 A THREE INTEGRAL DISTRIBUTION 
FUNCTION 

In this Section, we will construct a three-integral model us- 
ing a quadratic programming technique. With the three in- 
tegral modeling, we wish to approach two questions 

(i) Is it possible to fit the data without a central dark 
mass, with a distribution function that has the appropri- 
ate three integral structure? In this context, we recall the 
well-known case of M87, where Binney & Mamon (1982) 
succeeded in fitting the (then available) data with a distri- 
bution function that was very anisotropic. 

(ii) What is the phase-space structure of the best global 
31 model? 

This is an opportunity to present the use of the Quadratic 
Programming technique on a complex object, using a com- 
plete set of two-dimensional kinematical data. It should be 
clear from the outset, however, that we chose to produce an- 
alytical distribution functions. These implies that our third 
integral must be approximate. In that sense, the three in- 
tegral models we present here must be seen as belonging to 
the class of models with approximate third integrals. 



6.1 The method 

The essence of the method relies on the fact that most, if not 
all, observable kinematic quantities can be expressed as lin- 
ear functionals of the distribution function, if the potential 
is given. This means that if we approximate the distribu- 
tion function / as a sum of some conveniently chosen basis 
functions fi\ 



f 



Cifi{h,h,h), 



(2) 



then the same relation holds between the (e.g.) projected 
mass density p p and the projected mass density of the com- 
ponents p p ,i 



Pp = 



E 



Cj Pp , i 



or a similar relation for the projected pressures: 



(3) 



(4) 



If we stick to the mass density p for simplicity, we can con- 
struct a x 2 variable with the observed values pi as follows: 



Cipi,l 



(•») 



with p it i the value of the mass density at the I— th point for 
component i, and the wi associated weights. Clearly this \ 2 
is quadratic in the coefficients Ci. Inclusion of other observ- 
ables is trivial. 

The minimisation of such a \ 2 wn l produce a set of 
coefficients, which will yield a distribution function (Q) that 
will not necessarily be positive everywhere. Therefore, one 
must include constraints 



> C i F i [(7l,72,/ 3 )m] > 0, 



m = 1, 



(6) 



expressing the positivity of the distribution function on a 
grid in phase space with points {I\,l2, Iz)m- This formulates 
a problem of quadratic programming, which is a well-known 
problem for which efficient algorithms exist. 

In the present case, we use the luminosity model as pre- 
sented in section j^. As discussed in Sect. ^, a constant M/L 
model will not work at large radii. This is indicated by the 
rotation curve, which is almost flat up to 200" (Capaccioli 
et al. 1993). We therefore construct a model for the dark 
mass that preserves the form of the luminous mass (say, for 
simplicity of the argument, ellipses with fixed axis ratios and 
semi-major axes a) and that has the radial dependence 



/Otot(o) = Pium(a)(l + Ca p ). 

We used p — 1.5, and the coefficient C is such that 



ptot(a) 
plum (a) 



= 30, a = 30kpc. 



(7) 



(8) 



This prescription produces a flat rotation curve. In practice, 
we adopt the same shape for the dark matter as for the 
luminous mass density, by specifying the relation (Q) for all 
the radial functions in the harmonic expansion of pi um - In 
Table ^, we present the dependence of the mass to light ratio 
of the three-integral model on the radius in the equatorial 
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Table 4. Mass to light ratio of the three-integral model versus the 
radius in arcsec or kpc. The contribution of the luminous mass is 
given in the right column 
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8.7 
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46 
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38 
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35 
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63.1 


33 



plane 
Sect. |6.2 



For the sake of comparison, we will also present in 
a model where the mass to light ratio was forced 



to stay constant (besides the central dark mass). 

The calculation of the potential is straightforward, 
when the harmonic expansion of the total mass density is 
given. In practice, we interpolate the potential on a grid. 
Since it is our goal to explore what three integral dynam- 
ical model contributes to the modeling of NGC 3115, we 
need an expression for the third integral I3. This is done by 
fitting the potential with a Stackel potential, and adopting 
the third integral from the Stackel case as our approximation 
of the third integral. Thus, our models belong to the class 
of 31 models that use an approximation for the third inte- 
gral. The procedure is outlined in more detail in Dejonghe 
et al. (1996). 

Finally, we include a central singularity. It is a softened 
1/r singularity (hence a Plummer sphere), with a softening 
length of the order of 50000 AU (~ 0.25 pc or ~ 5 mas, 
about forty times smaller than the FOS aperture). In the 
transition region where the gravitational force of the black 
hole is of comparable strength as the gravitational force of 
the core of the galaxy, there is no third integral and the 
Stackel approximation of the potential is probably invalid. 
However, as argued in Sect. ^, we adopt as our working 
hypothesis that the spheroidal central component is two- 
integral. Hence, it will suffice if we produce three-integral 
models that do no penetrate too deep into the centre, where 
the model thus remains two-integral. It is however important 
to note that since orbits are non-local, the three-integral 
nature of our models significantly extends the probed space 
of solutions. 

The data we use to build the \ 2 are photometric and 
kinematical. As for the photometry, we produce a series of 
data points taken from a (logarithmically spaced) grid of the 
deprojected luminous mass density as obtained in Sect. [| 
The kinematical data for this modeling part are twofold. 

• Outside 50" along the major-axis and 8" along the 
minor-axis, the effect of pixel integration and seeing con- 
volution is negligible, and we assumed that the gaussian 
velocity and dispersion are good approximations for the two 
first true moments. We thus selected a set of points from the 



above mentioned sources. It is impractical (yet) to include 
all data points, while some are, mildly inconsistent. 

• The rest of the galaxy (mainly the central part and the 
major-axis inside 50") was sampled using the estimations of 
the t rue m oments given by the TIGER data presented in 
Sect. 2.2.1 and the published kinematics of vdM+94. This 



of course required to include the instrumental setup in the 

x 2 . 

As for the components, we use two families. The first 
one is an extention of Fricke components, and has distribu- 
tion functions of the form 



pmnEQ Jo 



= {E-E o f{J-J f 



(9) 



wherever ±J > Jo and E > Eo. They are zero elsewhere. 
The exponent m and n are an integer, but p can be real. In 
general the parameter p controls the central concentration of 
the component, the parameter q controls the amount of an- 
gular momentum, while n is indicative of the dependence on 
the third integral. The cut-off binding energy Eo is needed 
to limit the extent of the components, for further enhanced 
flexibility. Also, since the potential is almost singular at the 
centre, large values of p (thus components that are concen- 
trated) produce impossibly high mass densities at the centre. 
Hence, we rather control the central concentration with the 
parameter Eo, and adopt low or even fractional values for 
p. Finally, the cut-off angular momentum Jo prohibits the 
orbits to penetrate too deep into the potential well, where 
the third integral may not exist. 

The second family is designed to produce very thin 21 
discs in the equatorial plane of the galaxy (Batsleer & De- 
jonghe 1995), and has members of the form 



(s Z0 (j)rj 2 "(E-s Z0 (j)y 



(10) 



wherever ±J > 0, E > S Zo (J). They are zero elsewhere. 
The parameters p, q and s can be real. The quantity S zo (J) 
is the maximum binding energy that a star with angular 
momentum J can have if its distance from the equatorial 
plane remains below zq. If 20 = 0, So(J) is the binding 
energy of a circular orbit with J. The parameters p and q 
have the same meaning as in the first family, while s controls 
the concentration of the component towards the equatorial 
plane. 



6.2 Results 

In Fig. [Hi] we compare two models on the major axis, within 
the central 8". The solid line corresponds to a model with a 
6.5 x 10 8 MQblack hole, the dashed line has no black hole. 
Clearly, even in the three-integral case, the latter model fails 
dismally, especially in the inner 2". This is to be expected of 
course: the large dispersion in the centre exceeds the escape 
velocity if no black hole is present. Any reshuffling of orbits 
in order to obtain a favourable viewing position is insuffi- 
cient in the absence of a black hole, because there is simply 
not enough kinetic energy present. Although new compo- 
nents not included in our library may help, this is a strong 
suggestion that even three-integral models without a black 
hole cannot fit the data. 

In fact, if we compactify integral space by considering 
the representation in (E, J^/E, I3E) then powers, such as 
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Figure 16. Data points from B+94 (triangles) and vdM+96 
(crosses), and the best fit (31 models) with a 6.5 X 10 s MQblack 
hole (solid line) and best fit for a model without black hole (dot- 
ted line). The models have been convolved to correspond to B+94 
data. 




12 3 4 5 



u(kpc) 

Figure 17. A comparison of the luminous mass density, as ob- 
tained by the MGE technique (solid lines), and the 31 fit (dotted 
or dashed lines). The upper left panel shows the log of the lu- 
minous mass density along the major axis, the upper right panel 
dito along the minor axis 

Fricke components, form a complete set, if one allows nega- 
tive coefficients. In practice, of course, one is limited by finite 
computing resources. But at least, all orbits are present, with 
relative weights that vary continuously, and thus not nec- 
essarily with arbitrarily flexible relative weights. However, 
this practical incompleteness is not fundamentally different 
from the one present in, say, the Schwarzschild method. In 
the latter case, phase space cannot be covered completely, 
also for practical reasons, but the orbits considered can have 
arbitrarily flexible relative weights, at least without the (nec- 
essary) smoothing. So both methods are incomplete in prac- 
tice, be it to varying degree, but the incompleteness certainly 
will manifest itself differently. 



(6 
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Figure 18. The TIGER maps of the first true velocity moments 
V (upper panel) and a (lower panel), compared to the three- 
integral model with A/bh = 6.5 X 10 s M0(contours: step of 
20 km.s -1 for V and 25 km.s _1 for or). 



Regarding our second question, we first compare in 
Fig. [n] the luminous mass density as obtained with the MGE 
technique (see Section ^|) with the best 31 model that has a 
black hole of 6.5 x 10 8 M Q (see also Fig. §). The fit is quite 
reasonable, certainly considering the dynamic range that is 
to be covered (see upper panels). 

In Fig. Q [ii^ and |5(J we next compare the input kine- 
matical data with our best model that has a black hole of 
6.5 x 10 8 MQand includes a dark halo. The fit is good, ex- 
cept perhaps at the centre along the cut y = 20" and the one 
x = 40" (Fig. |(]). In fact, the data (Illingworth & Schechter 
1982) presented in this plot are mildly inconsistent with 
other published data along the minor or major axis (e.g. 
vdM+94, see Fig. |l9| ). The model with a constant mass to 
light ratio {M/L v = 10, dotted lines in Fig. |l| and|(]) does 
however show significant discrepancies along the major-axis 
and the cut y = 20" , with an overestimate of the velocity in 
the inner 50". This shows, as already stated, that the mass 
to light ratio must significantly increase outwards, by nearly 
a factor of 2 between the inner parts (R < 50", M/L ~ 6) 
and the last measured point (R ~ 200", M/L v ~ 10). We 
therefore confirm what was suggested by C+93 via simple 
dynamical arguments, but here using general three-integral 
dynamical models. Although a stellar mixture with a global 
M/Lv ~ 10 is not excluded, it would have to be recon- 
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Figure 19. A comparison of the kinematics along the major- 
axis and the best fits with a 6.5 X 10 8 M0black hole: including 
a dark halo (solid lines) or with a constant M/Ly = 10 (dotted 
lines). Inside 40 arcseconds we compare the model to the data 
of vdM+94, and outside this radius to the data of Capaccioli et 
al. (1993). The models have been convolved accordingly. 



ciled with the nearly constant (but decreasing) colours in the 
outer parts of NGC 3115 (for R > 120", Strom et al. 1977) 
and the metal poor populations observed in the halo by El- 
son (1997). The major-axis velocity and dispersion profiles 
are also more consistent with flat curves than with falling 
ones (see Fig. |l^ and Fig. |2o[ ) : this suggests that the M/Ly 
is still significantly increasing outwards for R > 150". In the 
view of these two points, we favour the hypothesis of the 
presence of a dark halo to explain the observed kinematics 
of NGC 3115 at large radii. 

For efficiency reasons, we included values for the her- 
mite coefficients at only 2 radii . At 7", the model convolved 
at the resolution of vdM+94 data has hs — —0.0025, close to 
the observed hs — —0.05, and in fact consistent with other 
published h 3 at this radius (Bender et al. 1994, F96). At 26" 
we obtain hz — —0.15, a much better fit than the 21 result, 
now perfectly consistent with the observed /13 = —0.16. In 
fact all observed higher order moments are now reasonably 
fit by the three-integral model. 

Finally, we present in Fig. [2l] and Fig. ^2] the mean 
rotation and velocity dispersions of the model in a merid- 
ian. The inner disc is tangentially anisotropic with a$/a r 
reaching a value of ~ 0.6 in the inner 2". This confirms 
the finding obtained via the two-integral model in Sect. ^ 
that the inner disc seems to be close to maximum rotation. 
Between 3" and 15" in the equatorial plane, the model is 
nearly isotropic. The inner spheroid (R < 60") is close to 
two-integral and has a nearly constant cr^/ov ~ 0.65. Go- 
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Figure 20. A comparison of the data points of Illingworth & 
Schechter (1982) and the best fits with a 6.5 X 10 8 M Q black 
hole: with a dark halo (see text, solid lines) and with a constant 
M/Ly = 10 (dotted lines). For each cut, we present the velocity 
(lower panel) and dispersion (upper panel). Note that the discrep- 
ancies between the model with the dark halo (solid lines) and the 
data are mostly due to inconsistencies in the data itself. 



ing outwards, the distribution function becomes significantly 
three- integral and more radially anisotropic. Already at a ra- 
dius of 30" in the equatorial plane (R ~ 1.5 kpc) a z /a r ~ 1.1 
and tx^/oy ~ 0.8. The ratio tx^/oy reaches a local maximum 
of 0.92, at a radius of ~ 100", where the outer disc contri- 
bution disappears and the dark matter becomes dominant. 



7 CONCLUSIONS 

We have presented a detailed analysis of the kinematics of 
NGC 3115 using different modeling techniques. As far as 
possible, we have made use of the photometric and kinemat- 
ical data available to us, in order to build dynamical models 
following realistic light and mass distributions. For the first 
time, we have included two-dimensional spectroscopic data 
in order to constrain the models in the central few hundred 
parcsecs. 

Jeans equations were used to define reasonable ranges 
for the mass to light ratio and the central dark mass in 
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Figure 21. Some moments of the best fit 31 distribution function 
in a meridian. Upper panel: (v) and lower panel: cr^. 
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Figure 22. Some moments of the best fit 31 distribution function 
in a meridian. Lower panel: oy and upper panel: cr z /cF r . 

the frame of a two integral dynamical model. Fitting the 
first two velocity moments in the central 40" gives Mbh = 
6.5±3.5x 10 8 M and M/L v = 6.5±0.7. These estimations 
were then tested by retrieving the corresponding distribu- 
tion functions via the Hunter & Qian method (Sect. ^). We 
would like to emphasize the fact that a simple two-integral 
model with a constant mass to light ratio can fit the velocity 
and dispersion profiles of NGC 3115 inside ~ 45" surpris- 
ingly well (see Sect. ^). 

However, no two-integral model could fit the value of 

V and /13 simultaneously^] as observed by vdM+94 around 
25" along the major-axis, where the outer disc significantly 
contributes to the surface brightness. Models built using the 
quadratic programming technique confirmed that three in- 

^ Remember that hz depends on V . 



tegral components solve this problem (Sect. g). The dynam- 
ical structure of the outer part of NGC 3115 is therefore 
very probably three-integral. 

In Sect. [B|, we showed that there are no two integral ax- 
isymmetric models which can fit the kinematics of NGC 3115 
in the central region, without the addition of a central dark 
mass of about 10 9 Mq. The best fit, which makes use 
of FOS/HST (Kormendy et al. 1996) and Gauss-Hermite 
moments, was obtained for Mbh = 0.94 x 10 9 MQand 
M/Lv = 6.1. This estimation of the central dark mass is 
at the upper limit of the range derived from simple Jeans 
models. This shows the need of high spatial resolution kine- 
matics including the higher order moments in order to more 
accurately constrain the mass distribution. 

The quadratic programming technique confirmed the 
need of a central dark mass as even the three-integral mod- 
els without a black hole failed to fit the central rise in the 
velocity dispersion. In principle, one could object that, since 
a library of 31 distribution functions is always limited, no 
matter what method is considered for the modeling, this is 
no definite proof that a roughly constant M/L model can- 
not fit the observables. However such a model would have 
to populate a lot of stars on radial orbits in the very centre. 
This would produce a galaxy core that is extremely prone 
to radial orbit instabilities, which, given the dynamical time 
of the order of 8 x 10 6 years at a radius of 10" (~ 500 pc), 
would have relaxed into a more stable configuration a long 
time ago, or would have produced a bar, evidence of which 
is not (yet) present. 

Previous estimations of the central dark mass in 
NGC 3115 are surprisingly different from one paper to the 
other. Kormendy & Richstone (1992) derived a value of Mbh 
from 1 to 2 x 10 9 M using spherical models corrected for 
flattening. Then Kormendy et al. (1996), using the same 
formalism, confirmed the upper limit of ~ 2 x 10 9 MQby 
including the SIS and HST/FOS data. Both papers used a 
mass to light ratio of M/Ly ~ 4. Finally, Magorrian et al. 
(1998) found a best fit value of M BH ~ 4.8 x 10 s M with 
a M/Lv ~ 6.75 (both values scaled for D — 10 Mpc), 
noting that their value for Mbh was probably an under- 
estimate. This is consistent with the picture drawn in the 
present work since we find that a value of Mbh as low 
as 6.5 x 10 s MQ(with a corresponding M/Lv = 6.8) is 
compatible with the existing data. However, our best fit 
model has Mbh = 0.94 x 10 9 Mq, and values as high as 
M BH = 2 x 10 9 M are excluded at more than the 3<r level. 
Our models, because they more precisely follow the observed 
light distribution and include the derivation of the full dis- 
tribution function, are improved estimations of the central 
dark mass and mass to light ratio. 

We finally showed that the mass to light ratio should in- 
crease by a factor of two between the inner parts (R < 50") 
and the outer parts (7? > 100") of NGC 3115. We produced 
a model which fits the flat rotation and dispersion profiles 
along the major-axis and thus gave the first evaluation of 
the amount of dark matter required to explain the observed 
kinematics at large radii: in our model the dark matter rep- 
resents about 50% of the total mass inside 150" and nearly 
70% inside 300". These values are obviously still uncertain 
as we did not fully examine alternative spatial distributions 
or radial profiles. 
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APPENDIX A: CONVOLUTION AND PIXEL 
BINNING USING A KERNEL 

It is possible to combine pixel integration and PSF smearing 
in only two integrals (as otherwise would be four: two for the 
seeing and two for the pixel integration). This is described 
in the following in the case of circular or rectangular pixels 
(see also Qian et al. 1995) 



Al Circular pixel 

For a circular pixel (TIGER like lenses) of radius d, the 
observable S at a position (x',y') after convolution with a 
single gaussian of dispersion a and integration on the pixel 
is: 



S a (x',y') 



where: 



1 

Txd 2 



oo p2n 



JO 



x rdr d6K(r) 



S[x' + r cos (6 + 6 ), y' + r sin (0 + O )] (Al) 



d/a 



•(f) 



.sr. _i( s _ T ./ CT )2 



(A2) 



K(r) = sdsl, 
Jo 

Note that if the integrals are derived with fixed quadratures, 
the Kernel has only to be calculated once, and only depends 
on the radius r. 



A2 Rectangular pixel 

For a rectangular pixel of dimension 21 x 2w. The observ- 
able So at a position (x',y') after convolution with a single 
gaussian of dispersion a and integration on the pixel is: 

poo p2ty 

S„(x',y') = — x / rdr d9K(r,9) 
4lw Jo Jo 

S [x + r cos (9 + 8a),y' + r sin {9 + 9 )] (A3) 



where: 
K(r,9) 



erf 



I + r cos ( 



V2a 



erf 



-I + r cos ( 



erf 



w + r sin ( 
V2a 



erf 



V2a 
w + r sin 
V2a 



(A4) 



and note again that if the integrals are derived with fixed 
quadratures, the Kernel can be calculated only once on a 
fixed grid. 



APPENDIX B: CUSPS IN MGE 
PHOTOMETRIC MODELS 

As we wish to build photometric models which can repro- 
duce central power laws, it has been necessary to extend the 
MGE formalism (Emsellem et al. 1996) to include such cusp 
components. The adopted form for the central component is 
then: 



v(m c ) = I c ■ 



2-ai 



-7/2 



x exp 



2-ai 



(Bl) 



where m 2 c — R 2 + Z 2 /q 2 . This is the natural generalized 
form of the 3D gaussian function. 

The integrated luminosity of this component is I c ■ 
2-7T h/2a c ) g c r (1.5 — 7/2), and the gravitational potential 
contribution simply: 



$(R,Z) = 4TrGq c aiP c x 

• rll-l/2,£?(R 2 + Z 2 /(l-e 2 c T 2 )) 



(1 - e 2 T 2 ) 



1/2 



dT 



(B2) 
(B3) 



with e 2 c = 1 - ql and P c = M/C x I c 
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